<Mark Rose>
<Section 2>
<12/4/18>
from matplotlib import pyplot as plt
from scipy.io import wavfile
import scipy.fftpack
import scipy.signal
import numpy as np
import IPython
import imageio
plt.rcParams["figure.dpi"] = 300 # Fix plot quality.
plt.rcParams["figure.figsize"] = (12,3) # Change plot size / aspect (you may adjust this).
class SoundWave(object):
"""A class for working with digital audio signals."""
# Problem 1.1
def __init__(self, rate, samples):
"""Set the SoundWave class attributes.
Parameters:
rate (int): The sample rate of the sound.
samples ((n,) ndarray): NumPy array of samples.
"""
self.rate = rate #Save the rate as an attribute
self.samples = samples #Save the samples as an attribute
return
raise NotImplementedError("Problem 1.1 Incomplete")
# Problems 1.1 and 1.7
def plot(self, force = False):
"""Plot the graph of the sound wave (time versus amplitude)."""
time = len(self.samples)/self.rate #Calculate the time by dividing the length of the samples devided by the rate
if force == False: #Plot the sound by itself if force is false
plt.plot(np.linspace(0,time, len(self.samples)), self.samples) #Plot the values according to time and the samples
plt.xlabel("Time (Seconds)")
plt.ylabel("Samples")
plt.title("Sound File Graph")
plt.ylim((-32768, 32767)) #Set the y limits
else: #Create two subplots if force is true
dft = scipy.fftpack.fft(self.samples) #Calculate the dft
freq = np.arange(len(dft))/len(dft)*self.rate #Calculate the frequencies
fig, axes= plt.subplots(1,2)
axes[0].plot(np.linspace(0, time, len(self.samples)), self.samples) #graph the original sound by doing time by samples
axes[0].set_xlabel("Time (Seconds)")
axes[0].set_ylabel("Samples")
axes[0].set_ylim((-32768,32767)) #Set the y limits
axes[0].set_title('Sound File Graph')
axes[1].plot(freq[0:len(freq)//2], np.abs(dft[0:len(freq)//2])) #Make a second subplot with half of the dft magnitues and frequencies
axes[1].set_xlabel('Frequency (Hz)') #Set the axis labels
axes[1].set_ylabel('Magnitude')
axes[1].set_title('Graph of Discrete Fourier Transform')
return
raise NotImplementedError("Problem 1.1 Incomplete")
# Problem 1.2
def export(self, filename, force=False):
"""Generate a wav file from the sample rate and samples.
If the array of samples is not of type np.int16, scale it before exporting.
Parameters:
filename (str): The name of the wav file to export the sound to.
"""
scaled = self.samples
real = np.real(self.samples)
if self.samples.dtype != 'int16' or force == True: #If the samples are not of type int16 or if the boolean force is True, scale the values
scaled = np.int16((real*32767.0/max(abs(real))))
wavfile.write(filename, self.rate, scaled) #Write a new sound file with the given filename, same rate, and new scaled values
return
raise NotImplementedError("Problem 1.2 Incomplete")
# Problem 1.4
def __add__(self, other):
"""Combine the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to add
to the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the combined samples.
Raises:
ValueError: if the two sample arrays are not the same length.
"""
if len(self.samples) != len(other.samples): #If the lengths are not the same raise an error
raise ValueError("Cannot add because the sample arrays from the SoundWave objects are not the same size")
combo = self.samples + other.samples #Add the two sounds together
return SoundWave(self.rate, combo) #Return a new SoundWave object with the two sounds and the original rate
raise NotImplementedError("Problem 1.4 Incomplete")
# Problem 1.4
def __rshift__(self, other):
"""Concatentate the samples from two SoundWave objects.
Parameters:
other (SoundWave): An object containing the samples to concatenate
to the samples contained in this object.
Raises:
ValueError: if the two sample rates are not equal.
"""
if self.rate != other.rate: #If the rates are not equal, raise an error
raise ValueError("The two sample rates are not equal")
combo = np.hstack((self.samples, other.samples)) #Stack the two vectors on top of each other
return SoundWave(self.rate, combo) #Return a sound wave object with the original rate and new combination
raise NotImplementedError("Problem 1.4 Incomplete")
# Problem 2.1
def __mul__(self, other):
"""Convolve the samples from two SoundWave objects using circular convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
my_samples = self.samples
other_samples = other.samples
if self.rate != other.rate:
raise ValueError("The two sample rates are not equal") #If the rates are not equal, raise an error
if len(self.samples) > len(other.samples):
num = len(self.samples) - len(other.samples)
padding = np.zeros(num)
other_samples = np.hstack((other_samples, padding))
elif len(other.samples) > self.samples:
num = len(other.samples) - len(self.samples)
padding = np.zeros(num)
my_samples = np.hstack(my_samples, padding)
conv = scipy.fftpack.ifft(scipy.fftpack.fft(my_samples)*scipy.fftpack.fft(other_samples)).real
return(SoundWave(self.rate, conv))
raise NotImplementedError("Problem 2.1 Incomplete")
# Problem 2.2
def __pow__(self, other):
"""Convolve the samples from two SoundWave objects using linear convolution.
Parameters:
other (SoundWave): An object containing the samples to convolve
with the samples contained in this object.
Returns:
(SoundWave): A new SoundWave instance with the convolved samples.
Raises:
ValueError: if the two sample rates are not equal.
"""
if self.rate != other.rate:
raise ValueError("The two sample rates are not equal") #If the rates are not equal, raise an error
n = len(self.samples)
m = len(other.samples)
i=0
while 2**i < n+m-1:
i+=1
pad1 = np.zeros(2**i-n)
pad2 = np.zeros(2**i-m)
my_samples = np.hstack((self.samples, pad1))
other_samples = np.hstack((other.samples, pad2))
conv = scipy.fftpack.ifft(scipy.fftpack.fft(my_samples)*scipy.fftpack.fft(other_samples)).real
return(SoundWave(self.rate, conv[0:n+m-1]))
raise NotImplementedError("Problem 2.2 Incomplete")
# Problem 2.4
def clean(self, low_freq, high_freq):
"""Remove a range of frequencies from the samples using the DFT.
Parameters:
low_freq (float): Lower bound of the frequency range to zero out.
high_freq (float): Higher boound of the frequency range to zero out.
"""
k_low = int(low_freq*len(self.samples)/self.rate)
k_high = int(high_freq*len(self.samples)/self.rate)
n = len(self.samples)
dft = scipy.fftpack.fft(self.samples)
dft[k_low:k_high] =0
dft[n-k_high: n-k_low]=0
self.samples = scipy.fftpack.ifft(dft).real
return
raise NotImplementedError("Problem 2.4 Incomplete")
SoundWave.__init__().SoundWave.plot().scipy.io.wavefile.read() and the SoundWave class to plot tada.wav.rate, samples = wavfile.read('tada.wav') #Get the rate and samples from the file
tada = SoundWave(rate, samples) #Create a SoundWave object
tada.plot() #Plot the resulting SoundWave
SoundWave.export().export() method to create two new files containing the same sound as tada.wav: one without scaling, and one with scaling (use force=True).IPython.display.Audio() to embed the original and new versions of tada.wav in the notebook.IPython.display.Audio(filename="tada.wav") #Play the original tada sound
tada.export("tada_noforce.wav", False) #Create a new sound file
IPython.display.Audio(filename="tada_noforce.wav") #Play the new sound file
tada.export("tada_withforce.wav", True) #Amplify the original tada sound and create a new file
IPython.display.Audio(filename="tada_withforce.wav") #Play the amplified sound
generate_note().generate_note() to create an A tone that lasts for two seconds. Embed it in the notebook.def generate_note(frequency, duration):
"""Generate an instance of the SoundWave class corresponding to
the desired soundwave. Uses sample rate of 44100 Hz.
Parameters:
frequency (float): The frequency of the desired sound.
duration (float): The length of the desired sound in seconds.
Returns:
sound (SoundWave): An instance of the SoundWave class.
"""
sample_space = np.linspace(0,duration, 44100*duration) #Create a group of values from 0 to the duration
samples = np.sin(2*np.pi*sample_space*frequency) #Use the function to get those values
rate = 44100 #Use the predefined rate and samples to make a SoundWave object
note = SoundWave(rate, samples)
return(note) #Now return the note
raise NotImplementedError("Problem 1.3 Incomplete")
my_note = generate_note(440, 2) #Create an A note for 2 seconds
IPython.display.Audio(rate=my_note.rate, data= my_note.samples) #Play the A note
SoundWave.__add__().SoundWave.__rshift__().a_note = generate_note(440,3) #Create an A, C, and E sounds that all last 3 seconds
c_note = generate_note(523.25,3)
e_note = generate_note(659.25,3)
a_minor_chord = a_note+c_note+e_note #Create an A minor chord
IPython.display.Audio(rate=a_minor_chord.rate, data = a_minor_chord.samples) #Play the chord
a_note = generate_note(440,1) #Make an a, c, and e note that are all one second long
c_note = generate_note(523.25,1)
e_note = generate_note(659.25,1)
a_minor_arpeggio = a_note >> c_note >> e_note #Concatenate the three sounds
IPython.display.Audio(rate=a_minor_arpeggio.rate, data = a_minor_arpeggio.samples) #Play the three sounds one after another
simple_dft() with the formula for $c_k$ given below.np.allclose() to check that simple_dft() and scipy.fftpack.fft() give the same result (after scaling).def simple_dft(samples):
"""Compute the DFT of an array of samples.
Parameters:
samples ((n,) ndarray): an array of samples.
Returns:
((n,) ndarray): The DFT of the given array.
"""
n = len(samples)
matrix_vals = np.reshape(np.arange(n), (n,1)) * np.reshape(np.arange(n), (1,n)) #Broadcast to make a matrix for the l and k values
dft = (samples * np.exp(-2*np.pi*1j*matrix_vals/n)).sum(axis=1)/n #Use the above equation to calculate the dft
return(dft)
raise NotImplementedError("Problem 1.5 Incomplete")
for n in range(6,20):
samples = np.random.random(n) #Make a bunch of random samples
my_dft = simple_dft(samples) #Check my dft vs their dft
given_dft = scipy.fftpack.fft(samples)
print(np.allclose(my_dft*n, given_dft)) #Print the results
simple_fft().simple_dft(), simple_fft(), and scipy.fftpack.fft().
Print the runtimes of each computation.np.allclose() to check that simple_fft() and scipy.fftpack.fft() give the same result (after scaling).def simple_fft(samples, threshold=1):
"""Compute the DFT using the FFT algorithm.
Parameters:
samples ((n,) ndarray): an array of samples.
threshold (int): when a subarray of samples has fewer
elements than this integer, use simple_dft() to
compute the DFT of that subarray.
Returns:
((n,) ndarray): The DFT of the given array.
"""
def split(g,N):
n = len(g)
if n <= N:
return(n*simple_dft(g)) #Use the predifined function with a small g
else:
even = split(g[::2],N)
odd = split(g[1::2],N) #Split into evens and odds
z = np.zeros(n)
z = np.exp(-2*np.pi*1j*np.arange(n)/n) #calculate the exponential parts using broadcasting
m = n // 2
return(np.hstack((even + z[:m] * odd, even + z[m:]*odd))) #return the concatenation of the two arrays
return(split(samples, threshold)/len(samples)) #Start the recursion
raise NotImplementedError("Problem 1.6 Incomplete")
powers_of_two = [1,2,4,8,16,32,64,128] #Create a list of powers of two
for n in powers_of_two:
samples = np.random.random(n) #Generate random sample vectors
my_dft = simple_fft(samples)
given_dft = scipy.fftpack.fft(samples) #Calculate the samples using my function, and the predifined one
print(np.allclose(my_dft*n, given_dft)) #Check to see if the functions have the same result
samples = np.random.random(8192)
%time simple_dft(samples) #This line times the simple dft algorithm
%time simple_fft(samples) #This line times the simple fft algorithm
%time scipy.fftpack.fft(samples) #This line times the predifined dft algorithm
SoundWave.plot() so that it accepts a boolean. When the boolean is True, take the DFT of the stored samples and plot (in a new subplot) the frequencies present on the $x$-axis and the magnituds of those frequences on the $y$-axis. Only the display the first half of the plot, and adjust the $x$-axis so that it correctly shows the frequencies in Hertz.a_note = generate_note(440,2) #Create an A note
a_note.plot(True) #Plot that note
a_note = generate_note(440,2) #The graph on here and the one in the book
c_note = generate_note(523.25,2) #match when time = 2 seconds
e_note = generate_note(659.25,2)
a_minor_chord = a_note+c_note+e_note #Create a chord with the notes a, c, and e
a_minor_chord.plot(True) #Use the plotting method to plot the sounds
Use the DFT to determine the individual notes that are present in mystery_chord.wav.
mrate, msamples = wavfile.read('mystery_chord.wav') #Get the rate and the samples from the sound file
myst= SoundWave(mrate, msamples) #Create a SoundWave object
dft=scipy.fftpack.fft(msamples)
dft=np.abs(dft) #Get the dft and the magnitudes
freq = np.arange(len(dft))/len(dft)*mrate #calculate the frequencies
vals = np.argsort(dft[0:len(freq)//2])[-4:] #Take the largest 3 from the first half
sounds = freq[vals]
sounds #These are the frequencies of the sounds
The notes are D, C, G, and A
SoundWave.__mul__() for circular convolution.tada.wav.tada.wav and the white noise. Embed the result in the notebook.rate, samples = wavfile.read('tada.wav')
tada = SoundWave(rate, samples) #create a SoundWave object with the rate and samples of tada
white_noise = SoundWave(rate, np.random.randint(-32767,32767, rate*2, dtype=np.int16)) #Create some random white noise in a SoundWave object
my_noise = white_noise * tada #Convolute the two noises
print(my_noise.rate)
my_noise = my_noise >> my_noise
IPython.display.Audio(rate=my_noise.rate, data = my_noise.samples) #Give the sound
SoundWave.__pow__() for linear convolution.CGC.wav and GCG.wav using SoundWave.__pow__() and scipy.signal.fftconvolve().SoundWave.__pow__() and scipy.signal.fftconvolve() sound the same.cgc_rate, cgc_samples = wavfile.read('CGC.wav')
cgc = SoundWave(cgc_rate, cgc_samples)
gcg_rate, gcg_samples = wavfile.read('GCG.wav')
gcg = SoundWave(gcg_rate, gcg_samples)
%time new_sound = cgc**gcg #This line times the simple dft algorithm
%time defined_samples = scipy.signal.fftconvolve(cgc_samples, gcg_samples) #This line times the simple fft algorithm
IPython.display.Audio(rate=new_sound.rate, data = new_sound.samples)
IPython.display.Audio(rate=cgc_rate, data = defined_samples)
Use SoundWave.__pow__() or scipy.signal.fftconvolve() to compute the linear convolution of chopin.wav and balloon.wav.
Embed the two original sounds and their convolution in the notebook.
c_rate, c_samples = wavfile.read('chopin.wav')
chop = SoundWave(c_rate, c_samples)
b_rate, b_samples = wavfile.read('balloon.wav')
ball = SoundWave(b_rate, b_samples)
chopin_echo = chop ** ball
IPython.display.Audio(rate=c_rate, data = c_samples)
IPython.display.Audio(rate=b_rate, data = b_samples)
IPython.display.Audio(rate=chopin_echo.rate, data = chopin_echo.samples)
SoundWave.clean().noisy1.wav by filtering out frequencies from $1250$-$2600$ Hz. Embed the original and the cleaned versions in the notebook.noisy2.wav. Embed the original and the cleaned versions in the notebook.n1_rate, n1_samples = wavfile.read('noisy1.wav')
noisy1 = SoundWave(n1_rate, n1_samples)
n2_rate, n2_samples = wavfile.read('noisy2.wav')
noisy2 = SoundWave(n2_rate, n2_samples)
IPython.display.Audio(rate=noisy1.rate, data = noisy1.samples)
noisy1.clean(1250,2600)
IPython.display.Audio(rate=noisy1.rate, data = noisy1.samples)
IPython.display.Audio(rate=noisy2.rate, data = noisy2.samples)
noisy2.plot(force=True)
noisy2.clean(1300,4400) #I determined from the above graph to filter from 1300 to 4400
IPython.display.Audio(rate=noisy2.rate, data = noisy2.samples)
vuvuzela.wav by filtering bad frequencies out of the left and right channels individually.vuvu_rate, vuvu_samples = wavfile.read('vuvuzela.wav')
vuvu_left = SoundWave(vuvu_rate, vuvu_samples[:,0])
vuvu_right = SoundWave(vuvu_rate, vuvu_samples[:,1])
vuvu_left.clean(200,500)
vuvu_right.clean(200,500)
new_sample = np.vstack((vuvu_left.samples, vuvu_right.samples))
print("Original Signal")
IPython.display.Audio(rate=vuvu_rate, data =vuvu_samples.T )
print('New Signal')
IPython.display.Audio(rate=vuvu_rate, data =new_sample )
license_plate.png so that the year printed on the sticker in the bottom right corner of the plate is legible.im = imageio.imread('license_plate.png')
fft2 = scipy.fftpack.fft2(im)
plt.imshow(np.log(np.abs(fft2)),cmap='gray')
plt.show()
mask = np.log(np.abs(fft2)) >12
fft2[mask] = np.mean(fft2)
plt.imshow(np.log(np.abs(fft2)),cmap='gray')
plt.show()
new_im = scipy.fftpack.ifft2(fft2).real
plt.imshow(new_im, cmap='gray')
The year on the sticker is 2013